#include "stdafx.h"
#include "gpMath.h"
#include "gpError.h"


//-------------------------------------------------------------------
// RealSmall : Returns the smallest positive real
//-------------------------------------------------------------------
 double RealSmall()
{
	return DBL_MIN;
}

//-------------------------------------------------------------------
// Abs : Returns the absolute value of a real
//-------------------------------------------------------------------
 double RealAbs(const double Value)
{
	return fabs(Value);
}

//-------------------------------------------------------------------
// IsEqual : Returns true if two reals are equal
//-------------------------------------------------------------------
 bool RealIsEqual(const double Value1, const double Value2)
{
	return RealAbs((Value1 - Value2)) < RealSmall();
}

//-------------------------------------------------------------------
// IsSimilar : Returns true if two reals are equal
//-------------------------------------------------------------------
 bool RealIsSimilar(const double One, const double Two)
{
	return RealIsEqual(One, Two);
}

//-------------------------------------------------------------------
// RealDigit : Returns the number of digits of precision in a real
//-------------------------------------------------------------------
 int RealDigits()
{
	return DBL_DIG;
}

//-------------------------------------------------------------------
// RealEpsilon : Returns the minimum positive real such that 
//               1.0 + x is not equal to 1.0
//-------------------------------------------------------------------
 double RealEpsilon()
{
	return DBL_EPSILON;
}

//-------------------------------------------------------------------
// RealFirst : Returns the minimum negative value of a real
//-------------------------------------------------------------------
 double RealFirst()
{
	return -DBL_MAX;
}

//-------------------------------------------------------------------
// RealFirst10Exp : Returns the minimum value of exponent(base 10) of
//                  a real.
//-------------------------------------------------------------------
 int RealFirst10Exp()
{
	return DBL_MIN_10_EXP;
}

//-------------------------------------------------------------------
// RealLast : Returns the maximum value of a real
//-------------------------------------------------------------------
 double RealLast()
{
	return  DBL_MAX;
}

//-------------------------------------------------------------------
// RealLast10Exp : Returns the maximum value of exponent(base 10) of
//                 a real.
//-------------------------------------------------------------------
 int RealLast10Exp()
{
	return  DBL_MAX_10_EXP;
}

//-------------------------------------------------------------------
// RealMantissa : Returns the size in bits of the matissa part of a 
//                real.
//-------------------------------------------------------------------
 int RealMantissa()
{
	return  DBL_MANT_DIG;
}

//-------------------------------------------------------------------
// RealRadix : Returns the radix of exponent representation
//-------------------------------------------------------------------
 int RealRadix()
{
	return  FLT_RADIX;
}

//-------------------------------------------------------------------
// RealSize : Returns the size in bits of an integer
//-------------------------------------------------------------------
 int RealSize()
{
	return BITS(double);
}

//-------------------------------------------------------------------
// IntToReal : Converts an integer in a real
//-------------------------------------------------------------------
 double IntToReal(const int Value)
{
	return Value;
}

//-------------------------------------------------------------------
// ATan : Returns the value of the arc tangent of a real
//-------------------------------------------------------------------
 double ATan(const double Value)
{
	return atan(Value);
}

//-------------------------------------------------------------------
// Ceiling : Returns the smallest integer not less than a real
//-------------------------------------------------------------------
 double RealCeiling(const double Value)
{
	return ceil(Value);
}

//-------------------------------------------------------------------
// Cos : Returns the cosine of a real
//-------------------------------------------------------------------
 double Cos(const double Value)
{
	return cos(Value);
}

//-------------------------------------------------------------------
// Cosh : Returns the hyperbolic cosine of a real
//-------------------------------------------------------------------
 double Cosh(const double Value)
{
	return cosh(Value);
}

//-------------------------------------------------------------------
// Epsilon : The function returns absolute value of difference
//           between 'Value' and other nearest value of
//           double type.
//           Nearest value is choseen in direction of infinity
//           the same sign as 'Value'.
//           If 'Value' is 0 then returns minimal positive value
//           of double type.
//-------------------------------------------------------------------
 double Epsilon(const double Value)
{
	double aEpsilon;

	if (Value >= 0.0) {
		aEpsilon = NextAfter(Value, RealLast()) - Value;
	}
	else {
		aEpsilon = Value - NextAfter(Value, RealFirst());
	}
	return aEpsilon;
}

//-------------------------------------------------------------------
// Exp : Returns the exponential function of a real
//-------------------------------------------------------------------
 double Exp(const double Value)
{
	return exp(Value);
}

//-------------------------------------------------------------------
// Floor : Return the largest integer not greater than a real
//-------------------------------------------------------------------
 double Floor(const double Value)
{
	return floor(Value);
}

//-------------------------------------------------------------------
// IntegerPart : Returns the integer part of a real
//-------------------------------------------------------------------
 double IntegerPart(const double Value)
{
	return ((Value > 0) ? floor(Value) : ceil(Value));
}

//-------------------------------------------------------------------
// Log10 : Returns the base-10 logarithm of a real 
//-------------------------------------------------------------------
 double Log10(const double Value)
{
	return log10(Value);
}

 //-------------------------------------------------------------------
// Pow : Returns a real to a given power
//-------------------------------------------------------------------
 double RealPow(const double Value, const double P)
{
	return pow(Value, P);
}

 //-------------------------------------------------------------------
 // RealPart : Returns the fractional part of a real.
 //-------------------------------------------------------------------
 double RealPart(const double Value)
 {
	 return fabs(IntegerPart(Value) - Value);
 }






//-------------------------------------------------------------------
// RealToInt : Returns the real converted to nearest valid integer.
//             If input value is out of valid range for integers,
//             minimal or maximal possible integer is returned.
//-------------------------------------------------------------------
 int RealToInt(const double Value)
{
	// Note that on WNT under MS VC++ 8.0 conversion of double value less 
	// than INT_MIN or greater than INT_MAX to integer will cause signal 
	// "Floating point multiple trap" (OCC17861)
	return Value < INT_MIN ? INT_MIN
		: Value > INT_MAX ? INT_MAX
		: (int)Value;
}

//-------------------------------------------------------------------
// Round : Returns the nearest integer of a real
//-------------------------------------------------------------------
 double RealRound(const double Value)
{
	return IntegerPart(Value + (Value > 0 ? 0.5 : -0.5));
}

//-------------------------------------------------------------------
// Sin : Returns the sine of a real
//-------------------------------------------------------------------
 double Sin(const double Value)
{
	return sin(Value);
}

 //-------------------------------------------------------------------
 // Sinh : Returns the hyperbolic sine of a real
 //-------------------------------------------------------------------
 double Sinh(const double Value)
 {
	 return sinh(Value);
 }

//-------------------------------------------------------------------
// ASinh : Returns the hyperbolic arc sine of a real
//-------------------------------------------------------------------
 double ASinh(const double Value)
{
	double  res;
	//  Modified by Sergey KHROMOV - Mon Nov 11 16:27:11 2002 Begin
	//  Correction of the formula to avoid numerical problems.
	//	res = log(X + sqrt(X * X + 1));
	if (Value > 0.)
		res = log(Value + sqrt(Value * Value + 1));
	else
		res = -log(sqrt(Value * Value + 1) - Value);
	//  Modified by Sergey KHROMOV - Mon Nov 11 16:27:13 2002 End
	return res;
}

 //-------------------------------------------------------------------
 // Square : Returns a real to the power 2
 //-------------------------------------------------------------------
 double RealSquare(const double Value)
 {
	 return Value * Value;
 }

//-------------------------------------------------------------------
// Tan : Returns the tangent of a real
//-------------------------------------------------------------------
 double Tan(const double Value)
{
	return tan(Value);
}

//-------------------------------------------------------------------
// Tanh : Returns the hyperbolic tangent of a real
//-------------------------------------------------------------------
 double Tanh(const double Value)
{
	return tanh(Value);
}














// ------------------------------------------------------------------
// Abs : Returns the absolute value of an Integer
// ------------------------------------------------------------------
 int  IntAbs(const int Value)
{
	return Value >= 0 ? Value : -Value;
}

// ------------------------------------------------------------------
// Hascode : Computes a hascoding value for a given Integer
// ------------------------------------------------------------------
 int  IntHashCode(const int theMe, const int theUpper)
{
	//return (Abs (theMe) % theUpper) + 1;
	return ((theMe & 0x7fffffff) % theUpper) + 1;
}

// ------------------------------------------------------------------
// IsEqual : Returns true if two integers are equal
// ------------------------------------------------------------------
 bool  IntIsEqual(const int theOne, const int theTwo)
{
	return theOne == theTwo;
}

// ------------------------------------------------------------------
// Hascode : Computes a hascoding value for a given unsigned integer
// ------------------------------------------------------------------
 int UIntHashCode(const __int32 theMe, const int theUpper)
{
	return ((theMe & 0x7fffffff) % theUpper) + 1;
}

// ------------------------------------------------------------------
// IsEqual : Returns true if two integers are equal
// ------------------------------------------------------------------
 bool UIntIsEqual(const __int32 theOne, const __int32 theTwo)
{
	return theOne == theTwo;
}

// ------------------------------------------------------------------
// IsSimilar : Returns true if two integers are equal
// ------------------------------------------------------------------
 bool UIntIsSimilar(const int One, const int Two)
{
	return One == Two;
}

// ------------------------------------------------------------------
// IsEven : Returns true if an integer is even
// ------------------------------------------------------------------
 bool IsEven(const int Value)
{
	return Value % 2 == 0;
}

// ------------------------------------------------------------------
// IsOdd : Returns true if an integer is odd
// ------------------------------------------------------------------
 bool IsOdd(const int Value)
{
	return Value % 2 == 1;
}

// ------------------------------------------------------------------
// Max : Returns the maximum integer between two integers
// ------------------------------------------------------------------
 int  IntMax(const int Val1, const int Val2)
{
	return Val1 >= Val2 ? Val1 : Val2;
}

// ------------------------------------------------------------------
// Min : Returns the minimum integer between two integers
// ------------------------------------------------------------------
 int  IntMin(const int Val1, const int Val2)
{
	return Val1 <= Val2 ? Val1 : Val2;
}

// ------------------------------------------------------------------
// Modulus : Returns the remainder of division between two integers
// ------------------------------------------------------------------
 int  IntModulus(const int Value, const int Divisor)
{
	return Value % Divisor;
}

// ------------------------------------------------------------------
// Square : Returns the square of an integer
// ------------------------------------------------------------------
 int IntSquare(const int Value)
{
	return Value * Value;
}

// ------------------------------------------------------------------
// IntegerFirst : Returns the minimum value of an integer
// ------------------------------------------------------------------
 int IntegerFirst()
{
	return INT_MIN;
}

// ------------------------------------------------------------------
// IntegerLast : Returns the maximum value of an integer
// ------------------------------------------------------------------
 int IntegerLast()
{
	return INT_MAX;
}

// ------------------------------------------------------------------
// IntegerSize : Returns the size in digits of an integer
// ------------------------------------------------------------------
 int IntegerSize()
{
	return BITS(int);
}


// ------------------------------------------------------------------
// Hascode : Computes a hascoding value for a given real
// ------------------------------------------------------------------
int HashCode(const double me, const int Upper)
{
	Standard_OutOfRange_Raise_if(Upper < 1, "Try to apply HashCode method with negative or null argument.");
	union
	{
		double R;
		int I[2];
	} U;
	//  U.R = Abs(me); // Treat me = -0.0 ADN 27/11/97
	U.R = me;
	return HashCode((U.I[0] ^ U.I[1]), Upper);
}

// ------------------------------------------------------------------
// ShallowCopy : Makes a copy of a real value
// ------------------------------------------------------------------
double ShallowCopy(const double me)
{
	return me;
}

//-------------------------------------------------------------------
// ACos : Returns the value of the arc cosine of a real
//-------------------------------------------------------------------
double ACos(const double Value)
{
	Standard_OutOfRange_Raise_if((Value < -1.) || (Value > 1.), "");
	return acos(Value);
}

//-------------------------------------------------------------------
// ACosApprox : Returns the approximate value of the arc cosine of a real.
//              The max error is about 1 degree near Value=0.
//-------------------------------------------------------------------

 double apx_for_ACosApprox(const double x)
{
	return  (-0.000007239283986332 +
		x * (2.000291665285952400 +
			x * (0.163910606547823220 +
				x * (0.047654245891495528 -
					x * (0.005516443930088506 +
						0.015098965761299077 * x))))) / sqrt(2 * x);
}

double ACosApprox(const double Value)
{
	double XX;
	if (Value < 0.) {
		XX = 1. + Value;
		if (XX < RealSmall())
			return 0.;
		return M_PI - apx_for_ACosApprox(XX);
	}
	XX = 1. - Value;
	if (XX < RealSmall())
		return 0.;
	return apx_for_ACosApprox(XX);

	// The code above is the same but includes 2 comparisons instead of 3
	//   double xn = 1.+Value;
	//   double xp = 1.-Value;
	//   if (xp < RealSmall() || xn < RealSmall())
	//     return 0.;
	//   if (Value < 0.)
	//     return M_PI - apx_for_ACosApprox (xn);
	//   return apx_for_ACosApprox (xp);
}

//-------------------------------------------------------------------
// ASin : Returns the value of the arc sine of a real
//-------------------------------------------------------------------
double ASin(const double Value)
{
	Standard_OutOfRange_Raise_if((Value < -1 || Value > 1), "");
	return asin(Value);
}

//-------------------------------------------------------------------
// ATan2 : Returns the arc tangent of a real divide by an another real
//-------------------------------------------------------------------
double ATan2(const double Value, const double Other)
{
	Standard_OutOfRange_Raise_if((Value == 0. && Other == 0.), "null value");
	return atan2(Value, Other);
}

//-------------------------------------------------------------------
// Sign : Returns |a| if B >= 0; -|a| if b < 0.
//             from x in the direction y
//-------------------------------------------------------------------
double Sign(const double a, const double b)
{
	//==== We use the function "nextafter()" fom library "math.h" ==============
	if (b >= 0.0) {
		return RealAbs(a);
	}
	else {
		return (-1.0 * RealAbs(a));
	}
}

//==========================================================================
//===== The special routines for "IEEE" and differents hardwares ===========
//==========================================================================
union RealMap {
	double real;
	unsigned int map[2];
};

//--------------------------------------------------------------------
// HardwareHighBitsOfDouble :  
//    Returns 1 if the low bits are at end.   (exemple: decmips and ALPHA )
//    Returns 0 if the low bits are at begin. (exemple: sun, sgi, ...)
//--------------------------------------------------------------------
static int HardwareHighBitsOfDouble()
{
	RealMap MaxDouble;
	MaxDouble.real = DBL_MAX;
	//=========================================================
	// reperesentation of the max double in IEEE is 
	//      "7fef ffff ffff ffff"   for the big indiens.
	//      "ffff ffff 7fef ffff"   for the littel indiens.
	//=========================================================

	if (MaxDouble.map[1] != 0xffffffff) {
		return 1;
	}
	else {
		return 0;
	}
}

//--------------------------------------------------------------------
// HardwareLowBitsOfDouble :  
//    Returns 0 if the low bits are at end.   (exemple: decmips )
//    Returns 1 if the low bits are at begin. (exemple: sun, sgi, ...)
//--------------------------------------------------------------------
static int HardwareLowBitsOfDouble()
{
	RealMap MaxDouble;
	MaxDouble.real = DBL_MAX;
	//=========================================================
	// reperesentation of the max double in IEEE is 
	//      "7fef ffff ffff ffff"   for the big indiens.
	//      "ffff ffff 7fef ffff"   for the littel indiens.
	//=========================================================

	if (MaxDouble.map[1] != 0xffffffff) {
		return 0;
	}
	else {
		return 1;
	}
}

static int HighBitsOfDouble = HardwareHighBitsOfDouble();
static int LowBitsOfDouble = HardwareLowBitsOfDouble();

double NextAfter(const double x, const double y)
{
	RealMap res;

	res.real = x;

	if (x == 0.0) {
		return DBL_MIN;
	}
	if (x == y) {
		//=========================================
		//   -oo__________0___________+oo
		//               x=y
		//  The direction is "Null", so there is nothing after 
		//=========================================

	}
	else if (((x < y) && (x >= 0.0)) || ((x > y) && (x < 0.0))) {
		//=========================================
		//   -oo__________0___________+oo
		//        y <- x     x -> y
		//
		//=========================================
		if (res.map[LowBitsOfDouble] == 0xffffffff) {
			res.map[LowBitsOfDouble] = 0;
			res.map[HighBitsOfDouble]++;
		}
		else {
			res.map[LowBitsOfDouble]++;
		}
	}
	else {
		//=========================================
		//   -oo__________0___________+oo
		//        x -> y     y <- x
		//
		//=========================================
		if (res.map[LowBitsOfDouble] == 0) {
			if (res.map[HighBitsOfDouble] == 0) {
				res.map[HighBitsOfDouble] = 0x80000000;
				res.map[LowBitsOfDouble] = 0x00000001;
			}
			else {
				res.map[LowBitsOfDouble] = 0xffffffff;
				res.map[HighBitsOfDouble]--;
			}
		}
		else {
			res.map[LowBitsOfDouble]--;
		}
	}
	return res.real;
}


//-------------------------------------------------------------------
// ATanh : Returns the value of the hyperbolic arc tangent of a real
//-------------------------------------------------------------------
double     ATanh(const double Value)
{
	Standard_OutOfRange_Raise_if(((Value <= -1.) || (Value >= 1.)), "Illegal agument in ATanh");
	double res;
	res = log((1 + Value) / (1 - Value)) / 2;
	return res;
}

//-------------------------------------------------------------------
// ACosh : Returns the hyperbolic Arc cosine of a real
//-------------------------------------------------------------------
double     ACosh(const double Value)
{
	Standard_OutOfRange_Raise_if((Value < 1.), "Illegal agument in ACosh");

	double  res;
	res = log(Value + sqrt(Value * Value - 1));
	return res;
}

//-------------------------------------------------------------------
// Log : Returns the naturaOPl logarithm of a real
//-------------------------------------------------------------------
double     Log(const double Value)
{

	Standard_NumericError_Raise_if(Value <= 0., "Illegal agument in Log");

	return log(Value);
}
//-------------------------------------------------------------------
// Sqrt : Returns the square root of a real
//-------------------------------------------------------------------
double     Sqrt(const double Value)
{
	Standard_OutOfRange_Raise_if((Value < 0.), "Illegal agument in Sqrt");
	return sqrt(Value);
}



int CharToInt(const char me)
{
	Standard_ConstructionError_Raise_if(!std::isdigit((unsigned char)me), "");

	char S[2];
	S[0] = me;
	S[1] = 0;
	return atoi(S);
}

// ------------------------------------------------------------------
// CharToInt : Converts a string in an integer value
// ------------------------------------------------------------------

int CharToInt(const char* me)
{
	const int Len = strlen(me);
	for (int I = 0; I < Len; I++) 
	{
		Standard_ConstructionError_Raise_if(!std::isdigit((unsigned char)me[I]),"");
	}
	return atoi(me);
}

// ------------------------------------------------------------------
// ShallowCopy : Copy of an integer
// ------------------------------------------------------------------

int ShallowCopy(const int me)
{
	return me;
}


// ------------------------------------------------------------------
// NextPrime : Compute the first prime number greater or equal than an integer
// ------------------------------------------------------------------

#define VALUESNBR 4

long NextPrime(const long me)
{

	struct svalue {
		int signiaib;
		int nbr;
	};

	struct svalue values[VALUESNBR];
	long ia;
	long maxia;
	long ib[4];
	int n[4];
	// int signiaib[4] = { -1 , +1 , +1 , -1 } ;
	int signiaib[4];
	signiaib[0] = -1;
	signiaib[1] = 1;
	signiaib[2] = 1;
	signiaib[3] = -1;
	long remain;

	int nbvalues;
	int loop;
	int nindd;
	long minn;
	long maxvn;
	long premret = 0;

	Standard_OutOfRange_Raise_if((me < 0 || me >
#if defined (__alpha) || defined(DECOSF1)
		127149130704178201
#else
		2147483647
#endif
		),		
		"Try to apply NextPrime method with negative, null or too large value.");
	if (me <= 7) {
		if (me <= 1)
			return 1;
		else if (me <= 2)
			return 2;
		else if (me <= 3)
			return 3;
		else if (me <= 5)
			return 5;
		else if (me <= 7)
			return 7;
	}

	minn = (me - 1) / 6;              // n minimum
	while (6 * minn + 1 < me) {
		minn += 1;
	}

	maxia = long(sqrt((double)me) / 6 + 1);

	maxvn = minn + VALUESNBR;

	nbvalues = 0;
	for (nindd = 0; nindd < VALUESNBR; nindd++) {
		if (6 * (nindd + minn) - 1 < me) {
			values[nindd].nbr = 1;
			values[nindd].signiaib = -1;
			nbvalues += 1;
		}
		else {
			values[nindd].nbr = 0;
			values[nindd].signiaib = 0;
		}
	}

	for (ia = 1; ia <= maxia; ia++) {
		if (nbvalues == VALUESNBR * 2) {
			break;
		}
		remain = -VALUESNBR;
		ib[0] = (minn + ia - remain) / (6 * ia - 1);
		n[0] = int(6 * ia*ib[0] - ia - ib[0] - minn);
		ib[1] = (minn - ia - remain) / (6 * ia - 1);
		n[1] = int(6 * ia*ib[1] + ia - ib[1] - minn);
		ib[2] = (minn + ia - remain) / (6 * ia + 1);
		n[2] = int(6 * ia*ib[2] - ia + ib[2] - minn);
		ib[3] = (minn - ia - remain) / (6 * ia + 1);
		n[3] = int(6 * ia*ib[3] + ia + ib[3] - minn);
		for (loop = 0; loop < 4; loop++) {
			if (n[loop] >= 0 && n[loop] < VALUESNBR) {
				if ((values[n[loop]].nbr == 0) ||
					(values[n[loop]].signiaib == signiaib[loop])) {
					values[n[loop]].signiaib = -signiaib[loop];
					values[n[loop]].nbr += 1;
					if (values[n[loop]].nbr <= 2)
						nbvalues += 1;
				}
			}
		}
	}
	for (nindd = 0; nindd < VALUESNBR; nindd++) {
		if (values[nindd].nbr == 0) {
			if (me <= 6 * (nindd + minn) - 1) {
				premret = 6 * (nindd + minn) - 1;
				break;
			}
			else if (me <= 6 * (nindd + minn) + 1) {
				premret = 6 * (nindd + minn) + 1;
				break;
			}
		}
		else if (values[nindd].nbr == 1) {
			if (values[nindd].signiaib > 0) {
				if (me <= 6 * (nindd + minn) - 1) {
					premret = 6 * (nindd + minn) - 1;
					break;
				}
			}
			else {
				if (me <= 6 * (nindd + minn) + 1) {
					premret = 6 * (nindd + minn) + 1;
					break;
				}
			}
		}
	}

	if (premret != 0) {
		return premret;
	}

	return NextPrime(6 * (maxvn - 1) + 2);

}
static const int Primes[] = {
	101,  1009,  2003,  3001,  4001,  5003,  6007,  7001,  8009,  9001,
  10007,        12007,        14009,        16007,        18013,
  20011,               23003,               26003,               29009,
					   33013,                      37003,
		 41011,                             46021,
		 51001,                                    57037,
									 65003,
  100019 };  // catch the biggest

const int NbPrimes = 26; // does not include the biggest

int NextPrimeForMap(const int N)
{
	int i;
	for (i = 0; i < NbPrimes; i++)
		if (Primes[i] > N) break;
	return Primes[i];
}
